# Title: De-stereotyping Public Performance Evaluation
# Yixin Liu & Chengxin Xu 2022
# International Public Management Journal
# IPMJ Replication file for Study 2

library(ggplot2) # plot
library(pwr) # power
library(ggpubr) # plot
library(effsize) # effect size
library(multiwayvcov) # cluster se
library(lmtest) # cluster se
library(egg) # ggplot(theme_article)
library(cowplot) # ggplot output
library(Rmisc) # plot stat table function
library(stargazer) # latex
library(xtable) # latex
library(broman) # power

####Function####
# figure title box functions
element_textbox <- function(...) {
  el <- element_text(...)
  class(el) <- c("element_textbox", class(el))
  el
}

element_grob.element_textbox <- function(element, ...) {
  text_grob <- NextMethod()
  rect_grob <- element_grob(calc_element("strip.background", theme_bw()))
  
  ggplot2:::absoluteGrob(
    grid::gList(
      element_grob(calc_element("strip.background", theme_bw())),
      text_grob
    ),
    height = grid::grobHeight(text_grob), 
    width = grid::unit(1, "npc")
  )
}

# percentage function
pcentFun <- function(x) {
  res <- na.omit(x)
  100 * (sum(res) / length(res))
}


# Setup working directory
WD = getwd()
setwd(WD)

####Data####
s2<-read.csv("Study2_Liu&Xu_IPMJ.csv")

####Data transformation####
# transform data from individual level to profile level
s2se<-subset(s2,rand!="JE")
s2je<-subset(s2,rand=="JE")
s2se$id<-c(1:668)
s2je$id<-c(1:334)

s2se$perf=s2se$se1
s2se$send=s2se$se2

#####JE####
s2JE<-data.frame("observation" = 1:668,
                 "IndID" = rep(c(1:334),each=2),
                 "Profile" = rep(c(1,2),334),
                 "perf" = NA, "send" = NA,
                 "mode" = "JE",
                 "SAT" = NA,
                 "student" = rep(c("White","Black"),334))

######performance####
###for the 1st profile
for (i in 1:334){
  resp<-as.numeric(s2je[which(s2je[,28]==i),6])
  perf<-which(s2JE$IndID==i & s2JE$Profile==1)
  s2JE$perf[perf]<-resp
}
###for the 2nd profile
for (i in 1:334){
  resp<-as.numeric(s2je[which(s2je[,28]==i),7])
  perf<-which(s2JE$IndID==i & s2JE$Profile==2)
  s2JE$perf[perf]<-resp
}

######school choice####
###for the 1st profile
for (i in 1:334){
  resp<-as.numeric(s2je[which(s2je[,28]==i),8])
  send<-which(s2JE$IndID==i & s2JE$Profile==1)
  s2JE$send[send]<-resp
}
###for the 2nd profile
for (i in 1:334){
  resp<-as.numeric(s2je[which(s2je[,28]==i),9])
  send<-which(s2JE$IndID==i & s2JE$Profile==2)
  s2JE$send[send]<-resp
}

######SAT####
# white
for (i in 1:334){
  resp<-as.numeric(s2je[which(s2je[,28]==i),20])
  SAT<-which(s2JE$IndID==i & s2JE$Profile==1)
  s2JE$SAT[SAT]<-resp
}
# black
for (i in 1:334){
  resp<-as.numeric(s2je[which(s2je[,28]==i),21])
  SAT<-which(s2JE$IndID==i & s2JE$Profile==2)
  s2JE$SAT[SAT]<-resp
}

######covariates####
# gender
for (i in 1:334){
  resp<-as.numeric(s2je[which(s2je[,28]==i),11])
  sex<-which(s2JE$IndID==i)
  s2JE$sex[sex]<-resp
}
# race
for (i in 1:334){
  resp<-s2je[which(s2je[,28]==i),12]
  race<-which(s2JE$IndID==i)
  s2JE$race[race]<-resp
}
# age
for (i in 1:334){
  resp<-as.numeric(s2je[which(s2je[,28]==i),13])
  age<-which(s2JE$IndID==i)
  s2JE$age[age]<-resp
}
# state
for (i in 1:334){
  resp<-s2je[which(s2je[,28]==i),14]
  state<-which(s2JE$IndID==i)
  s2JE$state[state]<-resp
}
# parent
for (i in 1:334){
  resp<-as.numeric(s2je[which(s2je[,28]==i),25])
  parent<-which(s2JE$IndID==i)
  s2JE$parent[parent]<-resp
}
# income
for (i in 1:334){
  resp<-as.numeric(s2je[which(s2je[,28]==i),16])
  income<-which(s2JE$IndID==i)
  s2JE$income[income]<-resp
}
# education
for (i in 1:334){
  resp<-as.numeric(s2je[which(s2je[,28]==i),17])
  education<-which(s2JE$IndID==i)
  s2JE$education[education]<-resp
}
# ideology
for (i in 1:334){
  resp<-as.numeric(s2je[which(s2je[,28]==i),18])
  ideology<-which(s2JE$IndID==i)
  s2JE$ideology[ideology]<-resp
}
# duration
for (i in 1:334){
  resp<-as.numeric(s2je[which(s2je[,28]==i),1])
  duration<-which(s2JE$IndID==i)
  s2JE$duration[duration]<-resp
}
# atpass
for (i in 1:334){
  resp<-as.numeric(s2je[which(s2je[,28]==i),26])
  atpass<-which(s2JE$IndID==i)
  s2JE$atpass[atpass]<-resp
}
# mcpass
for (i in 1:334){
  resp<-as.numeric(s2je[which(s2je[,28]==i),27])
  mcpass<-which(s2JE$IndID==i)
  s2JE$mcpass[mcpass]<-resp
}

# reorder observation
s2JE$observation<-s2JE$observation+668
s2JE$IndID<-s2JE$IndID+668

#####SE####
s2SE<-data.frame("observation" = 1:668,
                 "IndID" = 1:668,
                 "Profile" = 1,
                 "perf" = s2se$perf,
                 "send" = s2se$send,
                 "mode" = "SE",
                 "SAT" = s2se$RandomSAT,
                 "student" = s2se$rand,
                 "sex" = s2se$sex,
                 "race" = s2se$race,
                 "age" = s2se$age,
                 "state" = s2se$state,
                 "parent" = s2se$par,
                 "income" = s2se$income,
                 "education" = s2se$education,
                 "ideology" = s2se$ideology,
                 "duration" = s2se$duration,
                 "atpass" = s2se$atpass,
                 "mcpass" = s2se$mcpass)

#####Combine####
# combine SE JE datasets
SEJE<-rbind(s2SE,s2JE)
SEJE$student[SEJE$student=="SE1"]<-"White"
SEJE$student[SEJE$student=="SE2"]<-"Black"

####Regression####
SEJE$mode <-as.factor(SEJE$mode)
SEJE$mode = relevel(SEJE$mode, ref="SE")
SEJE$student <-as.factor(SEJE$student)
SEJE$student = relevel(SEJE$student, ref="White")

##performance
# no covariate
mod1<-lm(perf~mode+student+mode*student, data = SEJE)
summary(mod1)
vcov_1 <- cluster.vcov(mod1, SEJE$IndID)
v1<-coeftest(mod1, vcov_1)
v1

# with covariate
mod1c<-lm(perf~mode+student+mode*student
          +SAT+sex+race+age+state+parent+income+education+ideology, 
          data = SEJE)
summary(mod1c)
vcov_1c <- cluster.vcov(mod1c, SEJE$IndID)
v1c<-coeftest(mod1c, vcov_1c)
v1c

##school choice
# no covariate
mod2<-lm(send~mode+student+mode*student, data = SEJE)
summary(mod2)
vcov_2 <- cluster.vcov(mod2, SEJE$IndID)
v2<-coeftest(mod2, vcov_2)
v2

# with covariate
mod2c<-lm(send~mode+student+mode*student
          +SAT+sex+race+age+state+parent+income+education+ideology, 
          data = SEJE)
summary(mod2c)
vcov_2c <- cluster.vcov(mod2c, SEJE$IndID)
v2c<-coeftest(mod2c, vcov_2c)
v2c

####Table 2####
stargazer(mod1,mod1c,mod2,mod2c,
          style = "apsr", type="html", out="Table2.html",
          se=list(v1[,2],v1c[,2],v2[,2],v2c[,2]),
          p=list(v1[,4],v1c[,4],v2[,4],v2c[,4]),
          covariate.labels=c("JExBlack","Mode: JE",
                             "Students' major Race: Black"),
          column.labels=c("Percieve performance",
                          "Percieve performance (covariates)",
                          "School Choice",
                          "School Choice (covariates)"),
          omit = c("sex","race","age","state","parent","income",
                   "education","ideology","SAT"),
          dep.var.labels.include = FALSE,
          intercept.bottom=TRUE,
          order = c("modeJE:studentBlack","modeJE","studentBlack"),
          add.lines = list(c("Covariates", 
                             "No","Yes","No","Yes")),
          keep.stat=c("adj.rsq","n"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          align=TRUE,no.space=TRUE)

####Subgroup SE-JE####
# including SAT
SE<-subset(SEJE,SEJE$mode=="SE")
JE<-subset(SEJE,SEJE$mode=="JE")

# SE
mod3<-lm(perf~student+SAT, data=SE)
summary(mod3)

mod3c<-lm(perf~student+SAT
          +sex+race+age+state+parent+income+education+ideology, data=SE)
summary(mod3c)

mod4<-lm(send~student+SAT, data=SE)
summary(mod4)

mod4c<-lm(send~student+SAT
          +sex+race+age+state+parent+income+education+ideology, data=SE)
summary(mod4c)

# JE
mod5<-lm(perf~student+SAT, data=JE)
summary(mod5)
vcov_5 <- cluster.vcov(mod5, JE$IndID)
v5<-coeftest(mod5, vcov_5)
v5

mod5c<-lm(perf~student+SAT
          +sex+race+age+state+parent+income+education+ideology, data=JE)
summary(mod5c)
vcov_5c <- cluster.vcov(mod5c, JE$IndID)
v5c<-coeftest(mod5c, vcov_5c)
v5c

mod6<-lm(send~student+SAT, data=JE)
summary(mod6)
vcov_6 <- cluster.vcov(mod6, JE$IndID)
v6<-coeftest(mod6, vcov_6)
v6

mod6c<-lm(send~student+SAT
          +sex+race+age+state+parent+income+education+ideology, data=JE)
summary(mod6c)
vcov_6c <- cluster.vcov(mod6c, JE$IndID)
v6c<-coeftest(mod6c, vcov_6c)
v6c

####Table 3####
# performance
stargazer(mod3,mod3c,mod5,mod5c,
          style = "apsr", type="html",out="Table3.html",
          se=list(NULL,NULL,v5[,2],v5c[,2]),
          p=list(NULL,NULL,v5[,4],v5c[,4]),
          covariate.labels=c("Students' major Race: Black",
                             "Students' average SAT"),
          column.labels=c("Separate Evaluation",
                          "Joint Evaluation"),
          column.separate=c(2,2),
          omit = c("sex","race","age","state","parent","income",
                   "education","ideology"),
          dep.var.labels.include = FALSE,
          intercept.bottom=TRUE,
          order = c("studentBlack","SAT"),
          add.lines = list(c("Covariates", 
                             "No","Yes","No","Yes")),
          keep.stat=c("adj.rsq","n"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          align=TRUE,no.space=TRUE)

####Table 4####
# school choice
stargazer(mod4,mod4c,mod6,mod6c,
          style = "apsr", type="html",out="Table4.html",
          se=list(NULL,NULL,v6[,2],v6c[,2]),
          p=list(NULL,NULL,v6[,4],v6c[,4]),
          covariate.labels=c("Students' major Race: Black",
                             "Students' average SAT"),
          column.labels=c("Separate Evaluation",
                          "Joint Evaluation"),
          column.separate=c(2,2),
          omit = c("sex","race","age","state","parent","income",
                   "education","ideology"),
          dep.var.labels.include = FALSE,
          intercept.bottom=TRUE,
          order = c("studentBlack","SAT"),
          add.lines = list(c("Covariates", 
                             "No","Yes","No","Yes")),
          keep.stat=c("adj.rsq","n"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          align=TRUE,no.space=TRUE)


####Figure 4####
##histogram plots
# performance
SEJE.perf.na.omit<-subset(SEJE,!is.na(perf))
df1 <- summarySE(SEJE.perf.na.omit, measurevar="perf", groupvars=c("mode","student"))
df1$mode <- factor(df1$mode, levels = c("SE","JE"))
df1$student<-ifelse(df1$student=="White","White School","Black School")

p1<-ggplot(df1, aes(x=mode, y=perf, fill=student, color=student)) + 
  geom_bar(stat="identity", width = 0.7, position=position_dodge()) +
  geom_errorbar(aes(ymin=perf-ci, ymax=perf+ci), width=.1,
                position=position_dodge(0.7)) +
  labs(x=element_blank(), y=element_blank(), 
       title="DV: Perceived Performance") +
  scale_x_discrete(labels=c("Separate Evaluation","Joint Evaluation"))+
  theme_article()
p1<-p1+coord_cartesian(ylim=c(50,80))
p1<-p1+theme(legend.title=element_blank(),
             legend.text = element_text(size=12),
             axis.title.y=element_text(size=12),
             axis.text=element_text(size=10),
             legend.justification=c(0,1), legend.position=c(0,1),
             strip.background = element_blank(),
             strip.text.x = element_blank(),
             plot.title = element_textbox(
               hjust = 0.5, margin = margin(t = 5, b = 5)))+
  guides(colour = guide_legend(reverse=F), fill = guide_legend(reverse=F))
p1<-p1+scale_color_manual(values=c("black","black"))

p1

# school choice
SEJE.send.na.omit<-subset(SEJE,!is.na(send))
df2 <- summarySE(SEJE.send.na.omit, measurevar="send", groupvars=c("mode","student"))
df2$mode <- factor(df2$mode, levels = c("SE","JE"))
df2$student<-ifelse(df2$student=="White","White School","Black School")

p2<-ggplot(df2, aes(x=mode, y=send, fill=student, color=student)) + 
  geom_bar(stat="identity", width = 0.7, position=position_dodge()) +
  geom_errorbar(aes(ymin=send-ci, ymax=send+ci), width=.1,
                position=position_dodge(0.7)) +
  labs(x=element_blank(), y=element_blank(),
       title="DV: School Choice") +
  scale_x_discrete(labels=c("Separate Evaluation","Joint Evaluation"))+
  theme_article()
p2<-p2+coord_cartesian(ylim=c(50,80))
p2<-p2+theme(legend.title=element_blank(),
             legend.text = element_text(size=12),
             axis.title.y=element_text(size=12),
             axis.text=element_text(size=10),
             legend.justification=c(0,1), legend.position=c(0,1),
             plot.title = element_textbox(
               hjust = 0.5, margin = margin(t = 5, b = 5)))+
  guides(colour = guide_legend(reverse=F), fill = guide_legend(reverse=F))
p2<-p2+scale_color_manual(values=c("black","black"))

p2

f4.plot<-ggarrange(p1,p2+theme(axis.ticks.y = element_blank(),
                               axis.text.y = element_blank()), 
                   ncol = 2, nrow = 1)

f4.plot<-plot_grid(f4.plot)
f4.plot
#ggsave(file = "Figure4.jpg", width = 6, height = 3, dpi = 666)

####Figure 5####
##scatter plots
SEscatter1<-ggplot(SE, aes(x=SAT, y=perf, color=student)) + 
  geom_point(size=0.2)+
  geom_smooth(method=lm,aes(fill=student))+
  labs(title="Separate Evaluation",y="Perceived Performance")+
  scale_fill_manual(values=c("#00BFC4","#F8766D"),labels = c("White School", "Black School"))+
  scale_color_manual(values=c("#00BFC4","#F8766D"),labels = c("White School", "Black School"))+
  guides(colour = guide_legend(reverse=T), fill = guide_legend(reverse=T))+
  theme_bw()
SEscatter1<-SEscatter1+theme(legend.text = element_text(size=12),
                             axis.title.y=element_text(size=12),
                             axis.text=element_text(size=10),
                             legend.justification = c("right", "bottom"),
                             legend.background = element_blank(),
                             legend.title = element_blank(),
                             legend.position=c(0.95,0.05),
                             strip.background = element_blank(),
                             strip.text.x = element_blank(),
                             plot.title = element_textbox(
                               hjust = 0.5, margin = margin(t = 5, b = 5)))

JEscatter1<-ggplot(JE, aes(x=SAT, y=perf, color=student)) + 
  geom_point(size=0.2)+
  geom_smooth(method=lm,aes(fill=student))+
  labs(title="Joint Evaluation",y="Perceived Performance")+
  scale_fill_manual(values=c("#00BFC4","#F8766D"),labels = c("White School", "Black School"))+
  scale_color_manual(values=c("#00BFC4","#F8766D"),labels = c("White School", "Black School"))+
  guides(colour = guide_legend(reverse=T), fill = guide_legend(reverse=T))+
  theme_bw()
JEscatter1<-JEscatter1+theme(legend.text = element_text(size=12),
                             axis.title.y=element_text(size=12),
                             axis.text=element_text(size=10),
                             legend.justification = c("right", "bottom"),
                             legend.background = element_blank(),
                             legend.title = element_blank(),
                             legend.position=c(0.95,0.05),
                             strip.background = element_blank(),
                             strip.text.x = element_blank(),
                             plot.title = element_textbox(
                               hjust = 0.5, margin = margin(t = 5, b = 5)))

SEscatter2<-ggplot(SE, aes(x=SAT, y=send, color=student)) + 
  geom_point(size=0.2)+
  geom_smooth(method=lm,aes(fill=student))+
  labs(y="School Choice")+
  scale_fill_manual(values=c("#00BFC4","#F8766D"),labels = c("White School", "Black School"))+
  scale_color_manual(values=c("#00BFC4","#F8766D"),labels = c("White School", "Black School"))+
  guides(colour = guide_legend(reverse=T), fill = guide_legend(reverse=T))+
  theme_bw()
SEscatter2<-SEscatter2+theme(legend.text = element_text(size=12),
                             axis.title.y=element_text(size=12),
                             axis.text=element_text(size=10),
                             legend.justification = c("right", "bottom"),
                             legend.background = element_blank(),
                             legend.title = element_blank(),
                             legend.position=c(0.95,0.05))

JEscatter2<-ggplot(JE, aes(x=SAT, y=send, color=student)) + 
  geom_point(size=0.2)+
  geom_smooth(method=lm,aes(fill=student))+
  labs(y="School Choice")+
  scale_fill_manual(values=c("#00BFC4","#F8766D"),labels = c("White School", "Black School"))+
  scale_color_manual(values=c("#00BFC4","#F8766D"),labels = c("White School", "Black School"))+
  guides(colour = guide_legend(reverse=T), fill = guide_legend(reverse=T))+
  theme_bw()
JEscatter2<-JEscatter2+theme(legend.text = element_text(size=12),
                             axis.title.y=element_text(size=12),
                             axis.text=element_text(size=10),
                             legend.justification = c("right", "bottom"),
                             legend.background = element_blank(),
                             legend.title = element_blank(),
                             legend.position=c(0.95,0.05))

scatter.plot<-ggarrange(SEscatter1+theme(axis.ticks.x = element_blank(),
                                         axis.text.x = element_blank(),
                                         legend.position = "none",
                                         axis.title.x = element_blank()),
                        JEscatter1+theme(axis.ticks.y = element_blank(),
                                         axis.text.y = element_blank(),
                                         legend.position = "none",
                                         axis.title.y = element_blank(),
                                         axis.ticks.x = element_blank(),
                                         axis.text.x = element_blank(),
                                         axis.title.x = element_blank()), 
                        SEscatter2+theme(legend.position = "none"),
                        JEscatter2+theme(axis.ticks.y = element_blank(),
                                         axis.text.y = element_blank(),
                                         axis.title.y = element_blank()),
                        ncol = 2, nrow = 2)
scatter.plot<-plot_grid(scatter.plot)
scatter.plot
#ggsave(file = "Figure5.jpg", width = 6, height = 6, dpi = 666)

####Appendix B####
#####B2####
##Characteristics of sample

# gender dummy
s2$female<-ifelse(s2$sex==0,1,0)
# race
s2$white<-ifelse(s2$race=="white",1,0)
s2$black<-ifelse(s2$race=="black",1,0)
s2$his<-ifelse(s2$race=="hispanic",1,0)
s2$asian<-ifelse(s2$race=="asian",1,0)
s2$other<-ifelse(s2$race=="other",1,0)
# age
s2$young<-ifelse(s2$age<=29,1,0)
s2$midage<-ifelse(s2$age>29 & s2$age<=49,1,0)
s2$old<-ifelse(s2$age>49,1,0)
# income
s2$low<-ifelse(s2$income==1,1,0)
s2$med<-ifelse(s2$income>=2 & s2$income<5,1,0)
s2$high<-ifelse(s2$income>=5,1,0)
# education
s2$hed<-ifelse(s2$education>=5,1,0)
# ideology
s2$con<-ifelse(s2$ideology>=4,1,0)
s2$lib<-ifelse(s2$ideology<=2,1,0)
s2$mod<-ifelse(s2$ideology==3,1,0)

# subgroup single and pair
s2.se1<-subset(s2, rand=="SE1")
s2.se2<-subset(s2, rand=="SE2")
s2.je<-subset(s2, rand=="JE")

## ANOVA F test (randomization check)
# gender
female.aov <- aov(female ~ rand, data = s2)
summary(female.aov)
female<-summary(female.aov)[[1]][,5][1]
male.aov <- aov(sex ~ rand, data = s2)
summary(male.aov)
male<-summary(male.aov)[[1]][,5][1]
# race
wht.aov <- aov(white ~ rand, data = s2)
summary(wht.aov)
wht<-summary(wht.aov)[[1]][,5][1]
blk.aov <- aov(black ~ rand, data = s2)
summary(blk.aov)
blk<-summary(blk.aov)[[1]][,5][1]
his.aov <- aov(his ~ rand, data = s2)
summary(his.aov)
his<-summary(his.aov)[[1]][,5][1]
asian.aov <- aov(asian ~ rand, data = s2)
summary(asian.aov)
asian<-summary(asian.aov)[[1]][,5][1]
other.aov <- aov(other ~ rand, data = s2)
summary(other.aov)
other<-summary(other.aov)[[1]][,5][1]
# age
young.aov <- aov(young ~ rand, data = s2)
summary(young.aov)
young<-summary(young.aov)[[1]][,5][1]
mid.aov <- aov(midage ~ rand, data = s2)
summary(mid.aov)
midage<-summary(mid.aov)[[1]][,5][1]
old.aov <- aov(old ~ rand, data = s2)
summary(old.aov)
old<-summary(old.aov)[[1]][,5][1]
# income
low.aov <- aov(low ~ rand, data = s2)
summary(low.aov)
low<-summary(low.aov)[[1]][,5][1]
med.aov <- aov(med ~ rand, data = s2)
summary(med.aov)
med<-summary(med.aov)[[1]][,5][1]
high.aov <- aov(high ~ rand, data = s2)
summary(high.aov)
high<-summary(high.aov)[[1]][,5][1]
# education
hed.aov <- aov(hed ~ rand, data = s2)
summary(hed.aov)
hed<-summary(hed.aov)[[1]][,5][1]
# parent
parent.aov <- aov(par ~ rand, data = s2)
summary(parent.aov)
parent<-summary(parent.aov)[[1]][,5][1]
# ideology
lib.aov <- aov(lib ~ rand, data = s2)
summary(lib.aov)
lib<-summary(lib.aov)[[1]][,5][1]
con.aov <- aov(con ~ rand, data = s2)
summary(con.aov)
con<-summary(con.aov)[[1]][,5][1]
mod.aov <- aov(mod ~ rand, data = s2)
summary(mod.aov)
mod<-summary(mod.aov)[[1]][,5][1]

#####Table B1####
des_stat<-matrix(NA,18,9)
rownames(des_stat)<-c("Female","Male",
                      "White","Black","Hispanic","Asian","Other",
                      "Age: 18-29","Age: 30-49","Age: >= 50",
                      "Income: < $25k",
                      "Income: $25k to $75k",
                      "Income: >= $75k",
                      "College degree",
                      "Conservative","Liberal","Moderate",
                      "Parenthood")
colnames(des_stat)<-c("Frequency","%","Frequency","%",
                      "Frequency","%","Frequency","%",
                      "P-value")

tol.sum<-round(apply(s2[,c(28,11,29:43,25)],
                     2,sum,na.rm=T),digits=0)
se1.sum<-round(apply(s2.se1[,c(28,11,29:43,25)],
                     2,sum,na.rm=T),digits=0)
se2.sum<-round(apply(s2.se2[,c(28,11,29:43,25)],
                     2,sum,na.rm=T),digits=0)
je.sum<-round(apply(s2.je[,c(28,11,29:43,25)],
                    2,sum,na.rm=T),digits=0)

tol.p<-round(apply(s2[,c(28,11,29:43,25)],
                   2,pcentFun),digits=0)
se1.p<-round(apply(s2.se1[,c(28,11,29:43,25)],
                   2,pcentFun),digits=0)
se2.p<-round(apply(s2.se2[,c(28,11,29:43,25)],
                   2,pcentFun),digits=0)
je.p<-round(apply(s2.je[,c(28,11,29:43,25)],
                  2,pcentFun),digits=0)

des_stat[,1]<-tol.sum
des_stat[,3]<-se2.sum
des_stat[,5]<-se1.sum
des_stat[,7]<-je.sum
des_stat[,2]<-tol.p
des_stat[,4]<-se2.p
des_stat[,6]<-se1.p
des_stat[,8]<-je.p
des_stat[,9]<-round(c(female,male,wht,blk,his,asian,other,
                      young,midage,old,low,med,high,hed,
                      con,lib,mod,parent),2)

b1<-xtable(x=des_stat, digits = c(0,0,0,0,0,0,0,0,0,2))
#print.xtable(B1,type="html",file="b1.html")

#####B3####
# manipulation check

SEJEmcat<-subset(SEJE, mcpass==1 & atpass==1)
SEmcat<-subset(SE, mcpass==1 & atpass==1)
JEmcat<-subset(JE, mcpass==1 & atpass==1)

#####Figure B1####
##histogram plots
df1.mcat <- summarySE(SEJEmcat, measurevar="perf", groupvars=c("mode","student"))
df1.mcat$mode <- factor(df1.mcat$mode, levels = c("SE","JE"))
df1.mcat$student<-ifelse(df1.mcat$student=="White","White School","Black School")

p1.mcat<-ggplot(df1.mcat, aes(x=mode, y=perf, fill=student, color=student)) + 
  geom_bar(stat="identity", width = 0.7, position=position_dodge()) +
  geom_errorbar(aes(ymin=perf-ci, ymax=perf+ci), width=.1,
                position=position_dodge(0.7)) +
  labs(x=element_blank(), y=element_blank(), 
       title="DV: Perceived Performance") +
  scale_x_discrete(labels=c("Separate Evaluation","Joint Evaluation"))+
  theme_article()
p1.mcat<-p1.mcat+coord_cartesian(ylim=c(50,80))
p1.mcat<-p1.mcat+theme(legend.title=element_blank(),
                       legend.text = element_text(size=12),
                       axis.title.y=element_text(size=12),
                       axis.text=element_text(size=10),
                       legend.justification=c(0,1), legend.position=c(0,1),
                       strip.background = element_blank(),
                       strip.text.x = element_blank(),
                       plot.title = element_textbox(
                         hjust = 0.5, margin = margin(t = 5, b = 5)))+
  guides(colour = guide_legend(reverse=F), fill = guide_legend(reverse=F))
p1.mcat<-p1.mcat+scale_color_manual(values=c("black","black"))

p1.mcat

df2.mcat <- summarySE(SEJEmcat, measurevar="send", groupvars=c("mode","student"))
df2.mcat$mode <- factor(df2.mcat$mode, levels = c("SE","JE"))
df2.mcat$student<-ifelse(df2.mcat$student=="White","White School","Black School")

p2.mcat<-ggplot(df2.mcat, aes(x=mode, y=send, fill=student, color=student)) + 
  geom_bar(stat="identity", width = 0.7, position=position_dodge()) +
  geom_errorbar(aes(ymin=send-ci, ymax=send+ci), width=.1,
                position=position_dodge(0.7)) +
  labs(x=element_blank(), y=element_blank(),
       title="DV: School Choice") +
  scale_x_discrete(labels=c("Separate Evaluation","Joint Evaluation"))+
  theme_article()
p2.mcat<-p2.mcat+coord_cartesian(ylim=c(50,80))
p2.mcat<-p2.mcat+theme(legend.title=element_blank(),
                       legend.text = element_text(size=12),
                       axis.title.y=element_text(size=12),
                       axis.text=element_text(size=10),
                       legend.justification=c(0,1), legend.position=c(0,1),
                       plot.title = element_textbox(
                         hjust = 0.5, margin = margin(t = 5, b = 5)))+
  guides(colour = guide_legend(reverse=F), fill = guide_legend(reverse=F))
p2.mcat<-p2.mcat+scale_color_manual(values=c("black","black"))

p2.mcat

f1.mcat.plot<-ggarrange(p1.mcat,p2.mcat+theme(axis.ticks.y = element_blank(),
                                              axis.text.y = element_blank()), 
                        ncol = 2, nrow = 1)

f1.mcat.plot<-plot_grid(f1.mcat.plot)
f1.mcat.plot
#ggsave(file = "B1.jpg", width = 6, height = 3, dpi = 666)

#####Figure B2####
##scatter plots
SEscatter1mcat<-ggplot(SEmcat, aes(x=SAT, y=perf, color=student)) + 
  geom_point(size=0.2)+
  geom_smooth(method=lm,aes(fill=student))+
  labs(title="Separate Evaluation",y="Perceived Performance")+
  scale_fill_manual(values=c("#00BFC4","#F8766D"),labels = c("White School", "Black School"))+
  scale_color_manual(values=c("#00BFC4","#F8766D"),labels = c("White School", "Black School"))+
  guides(colour = guide_legend(reverse=T), fill = guide_legend(reverse=T))+
  theme_bw()
SEscatter1mcat<-SEscatter1mcat+theme(legend.text = element_text(size=12),
                                     axis.title.y=element_text(size=12),
                                     axis.text=element_text(size=10),
                                     legend.justification = c("right", "bottom"),
                                     legend.background = element_blank(),
                                     legend.title = element_blank(),
                                     legend.position=c(0.95,0.05),
                                     strip.background = element_blank(),
                                     strip.text.x = element_blank(),
                                     plot.title = element_textbox(
                                       hjust = 0.5, margin = margin(t = 5, b = 5)))

JEscatter1mcat<-ggplot(JEmcat, aes(x=SAT, y=perf, color=student)) + 
  geom_point(size=0.2)+
  geom_smooth(method=lm,aes(fill=student))+
  labs(title="Joint Evaluation",y="Perceived Performance")+
  scale_fill_manual(values=c("#00BFC4","#F8766D"),labels = c("White School", "Black School"))+
  scale_color_manual(values=c("#00BFC4","#F8766D"),labels = c("White School", "Black School"))+
  guides(colour = guide_legend(reverse=T), fill = guide_legend(reverse=T))+
  theme_bw()
JEscatter1mcat<-JEscatter1mcat+theme(legend.text = element_text(size=12),
                                     axis.title.y=element_text(size=12),
                                     axis.text=element_text(size=10),
                                     legend.justification = c("right", "bottom"),
                                     legend.background = element_blank(),
                                     legend.title = element_blank(),
                                     legend.position=c(0.95,0.05),
                                     strip.background = element_blank(),
                                     strip.text.x = element_blank(),
                                     plot.title = element_textbox(
                                       hjust = 0.5, margin = margin(t = 5, b = 5)))

SEscatter2mcat<-ggplot(SEmcat, aes(x=SAT, y=send, color=student)) + 
  geom_point(size=0.2)+
  geom_smooth(method=lm,aes(fill=student))+
  labs(y="School Choice")+
  scale_fill_manual(values=c("#00BFC4","#F8766D"),labels = c("White School", "Black School"))+
  scale_color_manual(values=c("#00BFC4","#F8766D"),labels = c("White School", "Black School"))+
  guides(colour = guide_legend(reverse=T), fill = guide_legend(reverse=T))+
  theme_bw()
SEscatter2mcat<-SEscatter2mcat+theme(legend.text = element_text(size=12),
                                     axis.title.y=element_text(size=12),
                                     axis.text=element_text(size=10),
                                     legend.justification = c("right", "bottom"),
                                     legend.background = element_blank(),
                                     legend.title = element_blank(),
                                     legend.position=c(0.95,0.05))

JEscatter2mcat<-ggplot(JEmcat, aes(x=SAT, y=send, color=student)) + 
  geom_point(size=0.2)+
  geom_smooth(method=lm,aes(fill=student))+
  labs(y="School Choice")+
  scale_fill_manual(values=c("#00BFC4","#F8766D"),labels = c("White School", "Black School"))+
  scale_color_manual(values=c("#00BFC4","#F8766D"),labels = c("White School", "Black School"))+
  guides(colour = guide_legend(reverse=T), fill = guide_legend(reverse=T))+
  theme_bw()
JEscatter2mcat<-JEscatter2mcat+theme(legend.text = element_text(size=12),
                                     axis.title.y=element_text(size=12),
                                     axis.text=element_text(size=10),
                                     legend.justification = c("right", "bottom"),
                                     legend.background = element_blank(),
                                     legend.title = element_blank(),
                                     legend.position=c(0.95,0.05))

scatter.mcat.plot<-ggarrange(SEscatter1mcat+theme(axis.ticks.x = element_blank(),
                                                  axis.text.x = element_blank(),
                                                  legend.position = "none",
                                                  axis.title.x = element_blank()),
                             JEscatter1mcat+theme(axis.ticks.y = element_blank(),
                                                  axis.text.y = element_blank(),
                                                  legend.position = "none",
                                                  axis.title.y = element_blank(),
                                                  axis.ticks.x = element_blank(),
                                                  axis.text.x = element_blank(),
                                                  axis.title.x = element_blank()), 
                             SEscatter2mcat+theme(legend.position = "none"),
                             JEscatter2mcat+theme(axis.ticks.y = element_blank(),
                                                  axis.text.y = element_blank(),
                                                  axis.title.y = element_blank()),
                             ncol = 2, nrow = 2)
scatter.mcat.plot<-plot_grid(scatter.mcat.plot)
scatter.mcat.plot
#ggsave(file = "B2.jpg", width = 6, height = 6, dpi = 666)

